GENETIC EXPONENTIALLY FITTED METHOD FOR SOLVING 
MULTI-DIMENSIONAL DRIFT-DIFFUSION EQUATIONS 
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Abstract. A general approach was proposed in this article to develop high-order exponentially 
fitted basis functions for finite element approximations of multi-dimensional drift-diffusion equations 
for modeling biomolecular electrodiffusion processes. Such methods are highly desirable for achieving 
numerical stability and efficiency. We found that by utilizing the one-one correspondence between 
continuous piecewise polynomial space of degree fc+1 and the divergence-free vector space of degree k, 
one can construct high-order 2-D exponentially fitted basis functions that are strictly interpolative 
at a selected node set but are discontinuous on edges in general, spanning nonconforming finite 
element spaces. First order convergence was proved for the methods constructed from divergence- 
free Raviart-Thomas space RT^ at two different node sets. 
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1. Introduction. We will propose in this article a general approach to construct 
exponentially fitted methods for numerically solving the drift-diffusion equation 

-V ■ (DVu + D/3uV(j>) = /, (1.1) 

where u is the density of charged particles, D is the diffusion coefficient, f3 is a 
problem-dependent constant assumed positive in this study, <fi is the given electrostatic 
potential field. The source function / models the generation or recombination of the 
particles, and is assumed here to be independent of u. The drift-diffusion equation 
is widely applied in semiconductor device modeling, where u{x) is the density of 
charge carriers [18], and in biomolecular simulations, where u{x) can be the density 
of ions, ligands, lipids, or other diffusive molecules [21, 22, 36]. Generalizations of 
Eq.(l.l) to model quantum effects, finite particle sizes, particle-particle correlations 
and interactions have been made recently [15, 23, 36], in many cases by replacing 
the electrostatic potential with an effective potential that models these non-electric 
interactions. The drift-diffusion equations arising in biomolecular simulations are in 
general multidimensional due to the intrinsic 3-D structure of macromolecules. For 
drift diffusion in bulk solution or through ion channels 3-D drift-diffusion equations 
are usually adopted so the charge density can be solved at sufficient temporal and 
spatial accuracy. For lateral transport of charged particles on surfaces, such as lateral 
motion of lipids, lipid clusters, or proteins on membrane surfaces, 2-D drift-diffusion 
will be very useful for the reduction of dimensionality [19]. Surface drift-diffusion 
equation rather than 2-D drift-diffusion equation has to be used for lateral diffusion 
on smoothly curved surfaces [3, 2, 36], for which the gradient and divergence operators 
in (1.1) have to be replaced by surface gradient and surface divergence operators, 
respectively. This surface drift-diffusion equation can be solved using surface finite 
element methods, which share many properties with 2-D finite element methods. It 
is known that numerical solutions of drift-diffusion equations may suffer instability 
in the advection-dominated regimes that usually appear near the strongly charged 
molecular surfaces such as ion channels, deoxyribonucleic acid (DNA), or membrane 
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surfaces with proteins bound or in the vicinity. Numerical solutions of drift-diffusion 
equation in these biomolecular systems necessitates stable and efficient numerical 
methods. 

Exponentially-fitted methods, known originally as the Scharfetter-Gummel method 
[32], is a class of finite element methods for solving Eq.(l.l) through an exponential 
approximation to the density function u(x) on individual mesh elements. This is made 
possible by the introduction of the Slotboom variable 



The density function u(x) can then be exponentially fitted through locally constant 
or linear approximations of the density current J = De~^'Vp. With exponentially 
fitted methods the gradient of the electric potential field can be incorporated to the 
basis functions, and this makes the methods intrinsically stable against strong drift 
when the electric field is strong. This desirable feature has motivated many studies 
of the exponentially fitted methods [1, 32, 10, 16, 30, 11, 26, 5]. 

Exponential fitting for 1-D drift-diffusion equation is obtained by exactly solving a 
two-point boundary value problem locally in a divergence- free space of Lagrange basis 
functions [32]. The resultant nodal basis functions for u{x) involve Bernoulli functions. 
In multi-dimensional space, however, it is generally difficult to find a divergence-free 
space where the basis functions are interpolative at nodes and continuous on edges, for 
the local problem does not admit an analytical solution on triangles, quadrilaterals, 
or tetrahedra. To circumvent this difficulty a number of averaging schemes, such as 
exact average, volume harmonic average, and edge harmonic average, were introduced 
to the exponentially transformed diffusion coefficient in the symmetrized equation 
(1.3). 2-D exponentially fitted methods of this type were constructed in [10]. Such 
a construction of multi-dimensional exponential fitting methods is further analyzed 
in [16], where it is found that the classical Scharfetter-Gummel exponentially fitted 
scheme can be exactly reproduced by the harmonic averaging of the piecewise linear 
Galerkin method. Linear convergence in the exponentially weighted norm is proved 
as well. While the convergence and stability of these multi-dimensional methods 
in numerical computations were not reported in all of these studies, it was shown 
that they can well resolve the sharp gradient of the charge density without incurring 
any unphysical oscillations. Averaging techniques, such as exact average and volume 
harmonic average, may suffer from overflow or lead to over-smeared solutions. Some of 
the best results are obtained using edge harmonic average, which utilizes approximated 
exponential fitting along edges of triangles. 

There are continual efforts to obtain genetic exponential fitting in multi-dimensional 
spaces rather than through averaging. These include approximate spline fitting [24, 
25, 30, 5, 4], nonconforming finite element [29], finite volume methods [20], and ex- 
ponentially fitted difference methods [17]. The approximate spline methods strongly 
resemble the 1-D exponential fitting as the locally divergence free current field solu- 
tions are sought. Such solutions, however, are not unique, and thus approximations 
will be in place to achieve uniqueness. Nevertheless, locally constant approximation 
to the density current in 2-D violates the nodal interpolation constraints [28]. This 
major drawback is overcome by using locally linear approximation to the density cur- 
rent, and nodal interpolative basis functions for density functions are obtained [31]. 



p = ue 



(1.2) 



with which Eq.(l.l) can be symmetrized to be 



(1.3) 
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These basis functions are discontinuous on edges in general so the finite clement space 
is nonconforming. Extension of this approach to 3-D equations seems feasible. The 
work in [5, 4] are among the few attempts to establishing 3-D exponential fitting 
methods, where 1-D exponential fitting is sought along individual edges connecting to 
the vertices, and the nodal basis functions are defined to be the linear combinations 
of local 1-D exponentially fitted basis functions. The resulting finite element space 
is conforming. Despite these efforts there are very few applications of exponentially 
fitted methods to realistic multi-dimensional drift-diffusion equations. The analysis 
and application in [26], for example, is on 1-D drift-diffusion equation with quantum- 
corrected potential. Sometimes Streamline-Upwind/Petrov-Galerkin (SUPG) meth- 
ods that were developed for advection-diffusion equations have to be adopted [12] for 
stabilization. Moreover, it seems unclear whether the current methodologies can be 
extended for constructing high order exponentially fitting methods. 

We are motivated to propose in this article a general approach for constructing ex- 
ponential fitting methods in multi-dimensional spaces. We start with the divergence- 
free vector space that is obtained by taking the curl of the nodal basis of k th order 
Raviart-Thomas conforming space Pk+i, to define a prototype of the basis functions 
for the Slotboom variable p. The final form of the exponentially fitted basis functions 
for p(x) and u{x) will be computed by enforcing the nodal interpolative constraints. 
Our approach is different from that in [31] as the latter seeks a divergence-free ap- 
proximation of the density current and restrictive assumptions are introduced to get 
analytical representations of the basis functions. We do not enforce the density cur- 
rent to be divergence-free, and this relaxation enables our approach to reproduce the 
standard Lagrange spaces for solving the Poisson equation at the limit of vanishing 
electric potential in the drift-diffusion equation. With our approach 2-D exponen- 
tially fitted methods of arbitrary high order can be readily constructed by starting 
with high-order Raviart-Thomas conforming spaces. We will show that our method 
docs not entail some of the crucial pitfalls indicated in previous studies. First-order 
3-D exponentially fitted method can be established similarly, but the extension of 
our approach to high-order methods for 3-D problems is difficult due to the lack of 
one-one correspondence between P^+i and the divergence-free vector space in 3-D. 

The rest of the paper is organized as follows. In the next section we review some of 
the most important exponential fitting methods and summarize their major features. 
In Section 2 we define the divergence- free spaces and their basis functions, and describe 
the procedure to compute exponentially fitted charge density u(x) using these basis 
functions. Here we shall show that our construction is general, can reproduce 2-D 
exponentially fitted method previously developed, and degenerates to the standard 
conforming Lagrange finite element spaces. In Section 3 first order convergence is 
proved for constructions based on RTq . We make conclusion remarks and outline the 
future work in Section 4. 

2. Genetic Multi-dimensional Exponential Fitting. We first introduce no- 
tations that will be used in the discussion below. Let fl € M 2 be a bounded domain 
with Lipschitz continuous boundary. We introduce for s > the standard Sobolev 
space H s (Tt) to the functions in Q, and use the standard inner product (-,-)s, the 
norm ||-|| s , and the semi-norm | • | s . The inner product and the norm in L 2 (fl) are 
denoted by || • | and (•,•)■ The subspace of L 2 (0) that consists of functions with zero 
mean value is denoted by Lq(TI). With these we define the space iJ(div; f2): 



#(div, n) := {v : v e (L 2 (ft)) 2 ; V • v e £ 2 (fi)} 
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equipped with the norm 

l|v|| ff(div;t2 ) = (||v|| 2 + ||V.v|| 2 ) 1/2 . 

On a regular triangulation Th of f2 we define the finite element space Vh 

V h = {v G H (div; Q) : v\ K € ^(K) 6 Tt; v ■ n| an = 0}, 

where n is the outer normal direction on the boundary dCl, and Vk is the space of 
vector- valued polynomials of degree k on the element K. To fulfill the condition of the 
vanishing normal component on the boundary dfl one usually chooses Vk to be the 
classical Raviart-Thomas element of order k (RTk) [27] or the Brezzi-Douglas-Marini 
element of order k (BDMk) [8]. For other iJ(div; f2) spaces constructed more recently 
the readers may refer to [35, 33]. 

We will work on the divergence- free subspace Dh of Vh, which is defined to be 

D h = {v G V h : (V • v, g) = V g G L 2 (Q)}. 

It is known [35, 9] that if Vk is chosen to be RTk then 

V ■ V h = {q e L 2 {n) : q\ K = P k {K)} 

or BDM k then 

V ■ V h = {q G £§(fi) : g|* = Pfc-i(iT)} 

where is the space of polynomials of degree < k, suggesting that Dh is indeed 
exactly divergence free: 

D h = {v G V h : V • v = 0}. 

In case that V/, = RTk, this subspace is also denoted RT^. 

According to the Helmholtz decomposition, any divergence-free vector v can be 
given as the curl of a potential function i/j: 

v = curh/> = (-ipy,tp x ): (2-1) 

where the subscripts denote the partial derivatives. Eq.(2.1) suggests that we can 
construct a divergence-free subspace by taking the curl of some appropriate space. 
This is actually possible, as stated by the following well-known result concerning the 
RTk and BDMk triangular elements [9, 34]: 

Theorem 2.1. There exists a one-to-one map curl: Sh — > Dh where the space 

S h = {V> e floHfi) : 1>k = Pk+x{K), K G T h } 

for triangular RTk elements with k > or BDMk elements with k > 1 . T/ie dimen- 
sion of Dh = RT® is equal to 

dimP k+1 (K) - 1 = i(A + l)(fc + 4). (2.2) 

Similar correspondence between iJ(curl) spaces and the divergence free subspaces 
Dh in 3-D is also indicated [6]. Two sample spaces of Dh = i?TjJ? on the reference 
triangle are given here: 

i?T o ° = span{(l,0) T ,(0,l) T }, (2.3) 
JZZ? - span{(x, 1 - 2a: - y) T , (0,-1 + Ax) T , {-x, yf , (1 - Ay, 0) T , (-1 + a; + 2y, -jfl^l) 
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We are now in a position to construct the exponentially fitted methods for solving 
the following boundary value problem for Eq.(1.3): Find p £ iiT 1 (fl) such that 

-V • J(p) = f in ft, 

J(p)=De- p +Vp = D(yu + 0uV<f>) 

p = ge M onr fl cffi, 

J(p)-n = on T N = dQ \ T D , 

where and Tn are the Dirichlet and Neumann subsets of the boundary dfl, re- 
spectively. We denote by Uh and Jh(uh) = D(\7iih + (3uh^7(f>) the finite element 
approximations of u and J(u). Let {v^} be the basis of on a triangular clement 
K, then the approximate density current function Jh{uh) in K can be given as a linear 
combination of {v^}: 

Jh{uh)\i< = ^2 QVi ' ( 2 - 6 ) 
»=i 

where is the dimension of Dh\x ■ Notice that Jh(uh) can be given in terms of the 
Slotboom variable ph = Uhe^^ as Jh(Ph) — De~^V ph, we shall have 

De- f) *Vp h \ K = Y,Wi- (2-7) 

»=i 

It remains to find the basis functions pj of the certain finite element space on K in 
which p/j can be approximated. 

2.1. Attempts to approximate density current exactly in divergence- 
free spaces. If each basis function pj has the representation 



Vft-n^^m^v,, (2.8) 



for some set of real numbers {rrij^}, then (i) Vp^ and Vph will be exactly divergence- 
free, and (ii) the right-hand side of Eq.(2.8), as the gradient of a scalar, is curl free. 
Consequently, its integration along any curve I (Pi, Pa) connecting two given points in 
K is path- independent, for 

\7 Pr dl = p J (P 2 )-p J (P 1 ). (2.9) 

i(Pi-Pi) 

For this reason one can choose the path to be the line connecting (x, y) and the starting 
point (xo,yo), c.f. Figure 2.1. Then the basis function pj(x,y) can be computed from 

i N k f S{x,y) 

Pj (x,y)=p j (x ,y ) + ^J2 rn ^ /„ e^^v^s), y(s)) ■ nds, (2.10) 



where S(x,y) — \/(x — xq) 2 + (y — yo) 2 and n = (x — xq, y — yo)/S is the directional 
vector pointing (xo,yo) to (x,y). Consequently, 

u j (x,y) = u j (x ,y ) + — ^m 3 ,, / e^^^\ t (x(s), y(s)) ■ nds, 



D 

i=l 







(2.H) 
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with Uj(xo,yo) = pj(xo,yo)e~^^ x °' v ■ These functions {uj(x,y)} form a basis func- 
tion of Uh in the element K . On the standard reference triangle one can choose a 
two-section path (xo,yo) — ► (x,yo) —> (x,y), each section parallel to the axes, giving 
rise to 



Pj (x,y) = Pj (x ,yo) + ^X>i,i ( C e ms ' Vo H{s,yo)ds+ f ' e?«*>%Kx,t)dt 

i = l \Jx Jyo 

(2.12) 

where (vf,v^) are the two components of the vector v,, and correspondingly 
Uj (x,y)=u j (x ,y )+-—J2 m ^ / e^°\f(s,y )ds+ / eP*^^{x,t)dJt) 

i—l \Jx Jyo J 

(2.13) 



3 




Fig. 2.1. One needs to choose the starting point (xq, yo) and the path from (xo,yo) to (x, y) 
for the integrals in Eq. (2. 10,2. 11). For RT® one can place (xo, yo) at one vertex and enforce 
interpolative constraints at all three vertices; or choose one of the edge centers be the start point 
and do enforcement at all three edge centers. Different basis functions will be generated for different 
paths. 



To uniquely determine the total Nk coefficients rriji and the constant Uj(xo,yo) 
in the expression (2.11) we will enforce the interpolative constraints at Nf. + 1 nodes. 
For example, for Dh = RTq on the triangle in Figure (2.1) we can expand equation 
(2.13) with j = 2 at three vertices (1, 2, 3) to get the following 



U2{x,y) = u 2 (x ,yo)-\ — / „ m M 



e ms < yo K*(s,yo)ds 



■I'n 



(2.14) 



If we choose (x±,yi) to be (xo,yo) then we must have from the interpolative 
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condition that 



p -/30(xi,j/i) N *> / rxi ryi \ 

= u 2 (ari,yi) = u 2 (xo, W )) + - Y, m *A / ^* ( ' ,so) vf (*, Vo)ds + / e^^vf {x l ,t\m>) 

U i—\ \Jx a Jyo J 

„-l3<t>(x 2 ,yi) Nk f r x z rV2 \ 

1 = u 2 (x 2 ,2/ 2 ) = «2(so,W>) + ^-E m 2,4 / eW s >y°K*(s,y )ds+ / e^^vf^^rt) 

j = i V^xo Jyo / 

= u 2 (x 3 ,y 3 )= u 2 (x ,y ) + ^m 2)j e^^vf ( S , 2/ )ds + / e^-'Vf (^3, ^fl7.) 

.-—1 \Jxq Jyn / 









1 




F12 


1 


F21 


F-22 


1 



These three equations form a linear system 

m 21 \ / 

m 22 = 1 I , (2.18) 

"2(2:0,2/0) / \ 

where Fu and i^i are defind as 

P -P<j>{x 2 ,y 2 ) r x 2 ryi 
F u = - n / e^ 8 ^vf(s,y )ds+ e M***\\ (x 3 , t)dt, 1 < 3(2.19) 

U J x Jyo 

p-/30(x3,j/ 3 ) rx 3 ry 3 

F 2l = - = / e^ s ^vf(s,y )ds + ^ x ^wV(x 3 ,t)dt, i< 3(2.20) 

Therefore, with the similar expansions of the other two basis functions U\(x, y), u 3 (x, y), 
one can arrive at the general linear system below: 

mn m 2 i rn 3 \ 

m 12 m 22 m 32 I = I, (2.21) 









1 




F12 


1 


F21 


F22 


1 



where I is the identity matrix, and 



e 



Fji = / e^^y^v^xis), y(s)) ■ nds, j < 2. (2.22) 

The following theorem proves that the matrix F is invertible for either set of the 
nodes, i.e., the vertices (1,2,3) or the edge centers (4,5,6). 

Theorem 2.2. The matrix F is nonsingular for the basis {v^} of RTq . 

Proof. We notice that every component Fji represents a linear transformation of 
Vj, i.e., Fji = Fj(\~i) where Tj : Dh\x ~> Consequently it is sufficient to prove that 
-^j( v ) = f° r J = 1, 2 if and only if v — 0. We only need to prove that if ^ (v) = 
for j = 1,2 then v = 0, as the reverse statement is trivial. We denote v by v = (a, 6) 
because any v 6 Dh\K is a constant vector. If the vertices are chosen to be the node 
set for enforcing the interpolative constraints, we have 

-P4>(x 2 ,y 2 ) rS(x 2 ,y 2 ) 

Ji(v) = / e^ x{s) ^ s) \an 2 - bn x )ds = p(an 2 - 6n x ), (2.23) 

D Jo 

-P<f>(x 3 ,y 3 ) fS(x 3 ,y 3 ) 

■Fa(v) = n / e^ ( - x ^' v ^(an 4 - bn 3 )ds = q(an 4 - bn 3 ), (2.24) 

D ./n 
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where (711,712) and (713,714) are the unit vectors in the direction (x\,yi) — > (£2,2/2) 
and (x\,yx) — > (£3,7/3), respectively, and p,q are constants. It following immediately 
that if Eqs. (2.23,2.24) are both equal to zero then a = b = since the vectors (n±, 712) 
and (773, 774) are not parallel. Similar result holds true if the edge nodes are chosen for 
enforcing the interpolative constraints, with (£, y) being the start point of the integral 
path. □ 

It turned out that parameters {rriji} determined as such do not fit the expansion 

N k 

(2.8). Indeed, if the curl- free nature of 771^ V; is enforced along with the interpola- 

i=l 

tive constraints we will have an over-determined problem for rriji , suggesting that one 
can not find a set of interpolative basis functions for Uh that constitute a constant 
density current. Since this curl-free condition is not satisfied, different basis func- 
tions Pj,Uj can be generated if different paths are chosen for the integrations. Figure 
(2.2) shows the difference between the two basis functions that are obtained along 
different integral paths. The differences on the edges suggest the basis functions on 
neighboring triangles are discontinuous in general. The finite element space is hence 
nonconforming. In addition, this method can not be generalized to construct high 
order exponentially fitted method. For example, direct computation on D] t = RT® in 
the reference triangle with <j> = 0, along the integration path (0, 0) — > (x, 0) — > (x, y), 
gives rise to a singular matrix F: 



( 1/8 





-1/8 


1/2 


-3/8 


1 \ 


1/2 





-1/2 


1 


-1/2 


1 


1/4 











-1/4 


1 


1/2 


-1 


1/2 





-1/2 


1 


3/8 


-1/2 


1/8 





-1/8 


1 


V 














1 / 




Fig. 2.2. Different basis functions Uj can be obtained by choosing different integral path s. Left: 
Difference of u\ in Eq.(2.11) computed along s being (0,0) — ¥ (x,0) —¥ (x,y) and (0,0) — > (0,y) — > 
(x,y). Right: Difference of u\ computed through (2.31) and (2.32). 

Remark 2.1. A divergence-free finite element space 

V)F(K) = span{l, e v ^ x , (V</> x x) • e 3 } 

is constructed in [ ], by assuming local linear approximations of the electric potential 
and the density current . The six parameters of the latter and one constant in Uj 
are determined by applying seven conditions, namely one divergence-free condition, 
three linear conditions arising in the curl-free nature ofVpu, and three interpolation 
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constraints. The linear approximation of the density current provides more degrees of 
freedom than the constant approximation in RT^ presented above. These additional 
degrees of freedom make it possible to enforce the curl-free condition along with the 
interpolative constraints. The component (V^xx)-e z of space allows the simulation of 
the diffusion flow orthogonal to the electric field lines. Interestingly, we can prove that 
the nonconforming space OF (if) is a special case of our construction. Consider the 
standard reference triangle and let D] x \k = RT® then = span{(l, 0) T , (0, 1) T }- 

Let Path 1 be defined by (1,0) —> (x,0) —> (x,y) and Path 2 be defined by (1,0) — > 
(0, y) —> (x, y). Note that by averaging the results for the two symmetric paths we will 
arrive at the final solution Uj(x,y) and if we assume the same linear potential in the 
element </>|x = ax + by + c, then for Path 1 we have 



~/34> N k / px py 

i—\ \Jx a Jyo 



Uj (x ,yo) 



D 

-P<t> 
~D~ 



•Hi 



(3a 



a f3(ax+c) 



(3a 



m i,2 



*_ e P(ax+by+c) ^_ e j3(ax+c) 



;3b 



(3b 



by) 



(3bD 



\(3aD (3bD ) PaD 



(2.25) 



For Path 2 we have 



Uj(x , yo) + ( m jA / e^ x °^ds + m j>2 / eP^dt 

'yo Jxa 



e~W 


E' 

i=l 


D 


e -W 


[m 


D 


e -P4> 


(m 


D 


m jt2 


+ e- 


(3aD 





*e w(lo ' s) <(0,s)((s- 

ya J xa 



Assuming that the integration in Eq. (2.22) is path-independent, then the results 
from Path 1 and Path 2 are the same, and equal to their average 



, , u]{x,y) + u 2 Jx,y) 
Uj{x,y) = - 1 ^ 

~ AX °' yo) + 2(3D\a + b) + \2f3aD 2(3bD ) 

+ e -fiax ( Jf^_ _ _ / _J_ p{ax+by) 

\2(3bD 2(3aD) jA \2f3bD 2[3aD ) V 1 
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Applying the approximation e x ~ 1 + x we arrive at 

Uj (», v )« ttj (^ lI/d ) + ^^- + -j + ^^-^j + ^-^j 

_ m .. (_L_ + _JL_ \ e -Ka X+by ) _ 0by _ JZ=vL V 8 ax - --^ 

J ' 1 l v 2/36Z? + 2/3a J D y ' POy \2paD 2pbD ) P<1X \2f3bD 2(3aD ) 

= a 3 + rjje-^-* + 7j (V0 x x) • e 3 , (2.28) 

for some reorganized constants atj,r]j,"/j. This last formula (2.28) is of the exact 
same form as in [ J. However, the restrictive assumption a, b =^ in the local linear 
representation of the electric potential is no longer necessary in our construction here 
because we do not compute the analytical integration of the function e' 3 * and therefore 
the analytical form of the potential <j) is not needed. When Uj of the representation 
in Eq.(2.27) is used to compute the current Jh, the first term will produce a current 
component in the direction of the electric field; the second term does not contribute 
to Jh, and the last two terms will provide additional fluxes that are not necessarily 
aligned with the electric field. Eq.(2.28) indicates that these additional fluxes are 
either not necessarily orthogonal to the electric field in general. Under very special 
circumstances, for instances a — b in Eq.(2.27), one can solve the same 777,^1,771^2, 
and this crosswind diffusion flux will disappear in Eq.(2.27) as it will be absorbed into 
the second term. 

It is worth noting the divergence-free condition was only approximately enforced 
to get a constant current field in [28]. 

2.2. General construction of exponentially fitted methods by giving up 
divergence-free condition. We would note that it is not fully justified to enforce the 
divergence- free condition in constructing the exponentially fitted finite element space. 
Divergence-free condition is seldom enforced in solving the Laplace equation, which 
is the limit of the drift-diffusion equation (1.1) at the limit of a vanishing electric 
field and with f(x) = 0. Divergence-free condition shall not be used for solving 
drift-diffusion equations with inhomogeneous source function f(x), or unsteady drift- 
diffusion equations. 

We will give up the divergence-free condition to seek an alternative construction. 
Our construction for the drift-diffusion equation must regenerate the standard finite 
element space Pk+i for solving the Poisson equation at the limit of vanishing electric 
field. These standard basis functions ipj of Pk+i can be derived from the known basis 
of Dh through the following correspondence 




(2.29) 



and the enforcement of the interpolative constraints. We are motivated to construct 
the basis pj for in a similar manner 




(2.30) 



On the reference triangle K, we can evaluate pj at an arbitrary point (x,y) via 
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different paths 
PA X 'V) = Pj( x o,Vo) 



-Y- 

D ^ 



flu 



5 M*.*)vf(a:, t)tft 

(2.31) 



1 Nk ( f 

Pj(x,y) = Pj(x ,y ) + Y)Y m 0' 1 (j 



(2.32) 

The expansion coefficients rriji and the constants pj(xo,yo) can be solved from the 
following linear system similar to (2.21) 



V o 



F 



1 



1 / 



V pi(xo,yo) 



m N k + l,l 



\ 



1 J 1 



m N k + lM h 

PN k (x ,y ) p Nk +i(x ,yo) ) 



l(N k + l)x(N k +l) 



(2.33) 

resulting from the enforcement of the interpolative constraints at selected node sets, 
with 



31 D 



in 



e m^,t) Y ^ x . :t ) dt \ f i < j < N k +i, i < i < N k , 



We can then evaluate 



Uj (x,y) = Pj (x,y)e-^y\ 



(2.34) 



(2.35) 



Proper scaling can make Uj(x,y) to be interpolative as well. It turns out in general 
that Nk + 1 = dimPfc +1 (_ftT), meaning the dimension of the finite element space for 
Uh\x corresponding to the divergence- free subspace D h \ K — RT^\ K is equal to the 
dimension of the space Su\k = Pk+i{K) whose curl produces Dh\K- The one-one 
correspondence between RT® and Pk+i indicates that the matrix F corresponding to 
Eq.(2.29), i.e., F with 



F — 



^ V i{s,yo)dt - 



Hi 



vf(xj,t)dt 



(2.36) 



must be invertible. We shall note that the finite element spaces constructed here are 
nonconforming and depend on the integral path, similar to the spaces obtained in 
subsection 2.1, because the solved parameters rriji may not fit Eq.(2.30). They do fit 
when = 0. 

Two examples of the node set for the finite element space u^Ik — spanjuj} 
constructed on Dh\ = RT§ on the reference triangle are shown in Figure 2.3. For a 

weak potential cf> = e^ 2 ^ x2+y2 the exponentially fitted basis functions are very similar 
to the basis of Pi(K) at the same set of nodes, c.f. Figure 2.4 and 2.6. The nature 
of exponential fitting is more clearly illustrated when the electric field is strong, c.f. 
Figure 2.5 and 2.7. It can be seen that the finite element spaces spanned by Pj, Uj are 
invariant with respect to affine linear maps, regardless of the choice of node set. This 
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feature is of particular importance to the estimation of interpolation errors through 
mappings to the reference triangles. The basis functions Uj constructed on — RT® 

with (f> = exp~ 4 ^/ x ' +y2 are shown in Figure 2.8, where the characters of interpolation 
at nodes and exponential fitting are also well depicted. The method developed here 
thus provides for the first time a feasible approach to systematically construct high- 
order exponentially fitted methods for solving 2-D drift-diffusion equation. 




Fig. 2.3. Difference node sets are chosen for enforcing the interpolative constraints for = 
RT® , giving rise to different basis functions of u^. Left: triangle vertices; Right: Edge nodes. The 
starting point (xq, j/q) and the path for the integration vary accordingly. 



1.2n 




Fig. 2.4. Basis functions Uj corresponding to RTq with interpolative constraints enforced at 
the three vertices. 4> = e~ 2 V x2 +y 2 . 

Our construction of the 2-D exponentially fitted methods can be illustrated by 
the left diagram in Figure 2.9. Such construction does not apply to 3-D problems 
because the mapping Pk+i ~> Qh is not one-to-one. Consequently, even a computable 
basis of Vh has been found and exponentially fitted, it will be difficult in general to 
compute the basis for Uh from these fitted basis functions, in contrast to 2-D cases. 
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0.8- \ 
0.6 




Fig. 2.5. Basis functions Uj corresponding to RT® with interpolative constraints enforced at 
the three vertices. <j> = 4e~ 2 ^/ x2 + y2 . 



1.5-, 




Fig. 2.6. Basis functions Uj corresponding to RT$ with interpolative constraints enforced at 
the three edge centers. <j> = e~ 2 V x2 +v 2 



It is possible however to follow the procedure indicated by (2.8) to construct 0(h) 
order 3-D exponentially fitted method, based on the fact that there is a simple basis 
of Dh\K = RTq I k in 3-D, that is v, = e, on a reference tetrahedron. Hence we can 



14 



MELISSA R. SWAGER AND Y. C. ZHOU 




use the expansion 

3 

Jh(uh)\K = y^QVi 

i=l 

and following the procedure identical to that in subsection 2.1 to finally get the 
following basis functions for Uh\^'- 

e -P4> 3 / r x rv 

u j (x,y,z)=u j (x ,y ,z )+---J2 m lA ^^^{s^z^ds + / e^ x ' t ^v v i {x,t,z Q )dt+ 

■ U i—i \Jx Jya 

I e^^ y ' r \t{x,y,r)dr J , (2.37) 

where v, = (vf , vf, vf ) T . The total four parameters, i.e., Uj(xo, yo, Zq) and the three 
coefficients m^j, will be determined by enforcing the interpolative constraints at four 
nodes, which can be the four nodes of the tetrahedral or the centers of its four faces. A 
linear system similar to Eq.(2.21) can be formed and its solvability can be established 
through a direct computation. The space span{uj} is nonconforming similar to the 
2-D cases. 

3. Convergence of Exponentially Fitted Methods based on BTq . In this 
section we shall give a convergence analysis of the Galerkin finite clement approxima- 
tion of the problem (2.5) using the nonconforming exponentially fitted basis functions 
(2.11) or (2.35). We consider the symmetrized form of the drift-diffusion equation, 
and adopt the approach in [31] to establish the 0(h) convergence of the exponentially 
fitted method constructed on Dh = RTq. Convergence of the high order methods 
and the 3-D method involve more complicated error estimate of the interpolation in 
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0.2 0.4 
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Fig. 2.8. Basis functions Uj corresponding to RT^ with interpolative constraints enforced at 
three vertices and three edge centers. <j> = 4e _2 v /a;2 +y 2 



curl 

p k+1 . +rt<; 



{«,}- 



integral 



B. 
a 

era 



Pk + i- 



erad curl 

- Qh ► v h 



integral integral 
Q;,— 



e^y 
— v h 



Fig. 2.9. Illustration of the construction of 2-D exponentially fitted methods (left). In 3-D 
(right), the mapping between Pk+i and is not one-to-one so constructing basis for from 
exponentially fitted Vh for arbitrary k is difficult. 



nonconforming finite element spaces and will be reported separately. Notice that the 
interpolative basis functions of the nonconforming finite element spaces for the Slot- 
boom variable ph can be obtained from Uj(x,y) through scaling. Let V = Hq(Q) the 
weak solution of the problem (2.5) can be uniquely solved from the following weak 
formulation 



(3.1) 
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where the coercive, symmetric bilinear form a(p, v) : V x V — > K is given 

a(p, v) = / J(p) ■ Vvdx = / De-^Vp ■ Vvdx. 
Jn Ja 

For the purpose of Galerkin finite element approximation of this weak formulation we 
define the finite element space 

V h = {we L 2 {9) : w e span{pj} V K € T; w = on <9Q}. 

By construction is continuous only at selected nodes but in general discontinuous on 
edges, hence Vh is not a subspace of V in general. Following the standard treatments 
for nonconforming finite element methods (E.g., Chapter 10 on [7]) we introduce the 
space V + Vh in which we can define discrete bilinear and linear form, as 

ah(ph,v h ) = V / Der^Vph ■ Vv h , V p h , v h eV + V h , (3.2) 
Jk 



Ken 



(f,Vh)h= Yl I f v *» VwfcGV},. (3.3) 
KeT h - " 



K 



We equip the space V + Vh with a norm 

1/2 



(3.4) 



It can be seen that this norm becomes the semi-norm |u|i,n for v £ V, which is itself 
a norm due to the Poincare inequality. To check that \\v\\h — if and only if v = for 
v € Vh we notice that v must be a constant for = 0. Since v € Vh is continuous 

at nodes by construction the constant must be the same for all elements, and hence 
shall be zero for v = on dfl. It follows immediately that a,h(phi v h) is continuous 
and coercive on V + Vh because 

\a h (v h ,w h )\= I De'^Vvh ■ Vw h < M\\v h \\h\\w h \\h, Vv h ,w h €V + V h , 
and 

a h (v h ,v h )=J2 [ De-^Vp h -Vv h >m\\v h \\l, Vv h eV + V h , 

where 

]\[ — e -/3min x( =n |0(x)| m _ p -f3 max xE n |0(x)| 

The nonconforming Galerkin finite element approximation of problem (1.3) finally 
reads: Find ph S Vh such that 

dh(ph,v h ) = (f,v h )h, V v h € Vh- (3.5) 



There are two major components in the convergence analysis of the nonconforming 
exponentially fitted finite element method, one is the interpolation error in Vh, and 
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the other is the consistency error arising from the nonconformity of the method, 
according to the second Strang lemma. The analysis in [31] follows this standard 
approach. Estimation of the interpolation error is obtained through a decomposition 

Hp - n hP \\ < ||p- n lP \\ h + \\n hP - n lP \\ h , (3.6) 

where TlhP is the interpolation of p in Vh, and Hip is the Pi conforming interpolation 
of p in Th- This decomposition is applicable to our constructions with interpolative 
constraints enforced at the same set of nodes as those for conforming Lagrange finite 
element spaces. If the edge nodes are chosen for the constraint enforcement we will 
have to adopt a decomposition alternative to (3.6) 

Hp - n,,p|| <\\ P - n lP \\ h + \\n hP - u lP \\ h , (3.7) 

here Hip is the P\ nonconforming interpolation of p in Th- The first terms in these 
two cases have a similar estimate ([13], P. 130 and [14]) 

||p — nip|| fc < Ci/»|p|a,n, ||p-nip|U<C2^p| 2) n (3.8) 

for different generic constants C\, C%- Consequently for both cases the following result 
concerning the convergence of the solution ph of Eq.(3.5) can be established. 

Theorem 3.1. Let p be the unique weak solution of Problem (1.3) in H 2 (Q) n 
Hq(£1) and let {Th} be a family of regular triangulations ofCl. Then there exist positive 
constants C\ (</>), C2(4>) and Cs{4>) that depend only on (f>, such that 

\\p-Ph\\h <C 1 (0)/ i |p| 2 , fi + (7 2 (^)||V^|| oo ,n|n 1 p| 1 ,n + C3^)/i||V0|| oo ,rJ|p||2 J n- (3.9) 
Proof. In [31]. □ 

4. Summary and Future Work. We proposed a general approach to construct 
exponentially fitted methods for solving 2-D drift-diffusion equations. Our construc- 
tions are based on the divergence-free subspace of -ff(div; Q), usually RT%. The basis 
functions Uj for the density function are computed from the basis of RT®, and through 
enforcement of interpolative constraints at selected node sets, giving rise to noncon- 
forming finite element spaces. We showed that our approach can reproduce some of 
the previous constructions of the exponentially fitted methods, and will recover the 
standard conforming or nonconforming finite element spaces when the electric poten- 
tial is zero. Thanks to the one-one mapping between RT® and Pk+i, our approach can 
be used to develop arbitrary high-order 2-D exponentially fitted methods from RT®. 
Extension of our approach for high-order 3-D methods is difficult because there does 
not exist a similar one-one correspondence between k th order divergence-free vector 
space and Pk+i, but construction of first-order 3-D exponentially fitted method is 
possible. 

General convergence theory for the high-order methods can be established fol- 
lowing decomposition (3.6), the estimate of Pf. conforming interpolation errors, the 
estimate of II^p — II^p, and the consistency error of high-order nonconforming finite 
element methods. These results will be reported in a forthcoming paper. We will also 
conduct extensive numerical experiments of the methods, in particular the high order 
methods, and apply them to realistic biomolecular drift-diffusion problems. 
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